library(cem)
library(tidyverse)

rm(list=ls())
load("~/Dropbox/PHD/Indonesia_GPS/data-for-analyses/virginclip_pop.RData")

virginclip_pop <- virginclip_pop %>%
  dplyr::select(-geometry) %>%
  filter(peat_depth_GW>0)

# outcome vars

virginclip_pop$did2016 <- virginclip_pop$forestha_2016 - virginclip_pop$forestha_2011
#virginclip_pop$did2005 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2005

virginclip_pop$did2017 <- virginclip_pop$forestha_2017 - virginclip_pop$forestha_2011
#virginclip_pop$did2004 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

virginclip_pop$did2018 <- virginclip_pop$forestha_2018 - virginclip_pop$forestha_2011
#virginclip_pop$did2003 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2003

t1 = 2010
t11 = 2011
t0 = 2004

t2 = 2016
t3 = 2017
t4 = 2018

virginclip_pop$didpre = virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

virginclip_pop$ddd2016 <- virginclip_pop$did2016 - ((t2-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2017 <- virginclip_pop$did2017 - ((t3-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2018 <- virginclip_pop$did2018 - ((t4-t11)/(t1-t0))*virginclip_pop$didpre

##### Coarsened Exact Matching #####

distroads_cut = seq(0,300,25)
agc_breaks3 = seq(0, max(virginclip_pop$agc_bcc_mean), 30)
agc_breaks4 = seq(0, max(virginclip_pop$agc_bcc_mean), 50)

##### Model 1 #####

vars1 = c("forestha_2005","forestha_2006", "forestha_2007", "forestha_2008", "forestha_2009", "forestha_2010",
           "slope", "elev", "pop_2005", "pop_2010", "distroads", "agc_bcc_mean", "dist_palms_km")

# imbalance(group=virginclip_pop$moratorium, data=virginclip_pop[vars])

palms_cut = seq(0, max(virginclip_pop$dist_palms_km), 50)

nvars1 = setdiff(colnames(virginclip_pop), vars1)
mat1 <- cem(treatment = "moratorium", data = virginclip_pop, drop = nvars1, 
             cutpoints = list(agc_bcc_mean = agc_breaks3, distroads = distroads_cut, 
                              dist_palms_km = palms_cut), keep.all=TRUE)

est1 = att(mat1, ddd2017 ~ moratorium, data = virginclip_pop)

##### Model 2 #####

vars2 = c("forestha_2005","forestha_2006", "forestha_2007", "forestha_2008", "forestha_2009", "forestha_2010",
           "slope", "elev", "pop_2005", "pop_2010", "distroads", "agc_bcc_mean")

# imbalance(group=virginclip_pop$moratorium, data=virginclip_pop[vars])

nvars2 = setdiff(colnames(virginclip_pop), vars2)
mat2 <- cem(treatment = "moratorium", data = virginclip_pop, drop = nvars2, 
             cutpoints = list(agc_bcc_mean = agc_breaks3, distroads = distroads_cut), keep.all=TRUE)

est2 = att(mat2, ddd2017 ~ moratorium, data = virginclip_pop)


##### Model 3 #####

vars3 = c("forestha_2005","forestha_2006", "forestha_2007", "forestha_2008", "forestha_2009", "forestha_2010",
           "slope", "elev", "pop_2005", "pop_2010", "distroads", "agc_bcc_mean", "dist_palms_km", "dist_timber_km",
           "dist_logging_km")

timbs_cut = seq(0, max(virginclip_pop$dist_timber_km), 50)
logs_cut = seq(0, max(virginclip_pop$dist_logging_km), 50)


# imbalance(group=virginclip_pop$moratorium, data=virginclip_pop[vars])

nvars3 = setdiff(colnames(virginclip_pop), vars3)
mat3 <- cem(treatment = "moratorium", data = virginclip_pop, drop = nvars3, 
             cutpoints = list(agc_bcc_mean = agc_breaks3, distroads = distroads_cut,
                              dist_palms_km = palms_cut, dist_timber_km = timbs_cut,
                              dist_logging_km = logs_cut), keep.all=TRUE)

est3 = att(mat3, ddd2017 ~ moratorium, data = virginclip_pop)

##### Model 4 #####

vars4 = c("forestha_2005","forestha_2006", "forestha_2007", "forestha_2008", "forestha_2009", "forestha_2010",
           "slope", "elev", "pop_2005", "pop_2010", "distroads", "agc_bcc_mean")

# imbalance(group=virginclip_pop$moratorium, data=virginclip_pop[vars])

nvars4 = setdiff(colnames(virginclip_pop), vars4)
mat4 <- cem(treatment = "moratorium", data = virginclip_pop, drop = nvars4, 
             cutpoints = list(agc_bcc_mean = agc_breaks4, distroads = distroads_cut), keep.all=TRUE)

est4 = att(mat4, ddd2017 ~ moratorium + dist_palms_km + dist_timber_km + dist_logging_km, data = virginclip_pop)



# remove data and save results
rm(virginclip_pop)
save.image("~/results/cem-results-peat.RData")
